快速傅里叶变换 FFT 笔记
引入
在学数论时会遇到多项式乘法。而快速傅里叶变换可以帮助你在 \(O(n log(n))\) 的时间复杂度内解决此问题，当然常数大不大另说。
前言
一个多项式大概长这样：
为了表示方便，我们用 \(F\) 表示这个多项式。
考虑将 \(F\) 看作一个函数，将一些值代入而获得 \(F\) 的一些点值。
这有什么用呢？
我们知道将一个函数的系数表达可以转化为点值表达，而从点值表达也可以转化为系数表达，后者具体方式形似解多元一次方程，只要有与未知数相同数量的点值，就得以解出系数表达式。
假设将要相乘的 \(F\) 、\(G\) 乘起来等于 \(H\) ，我们先求出 \(F\) 、 \(G\) 相应数量的点值，此时再将对应的点值相乘，得到 \(H\) 对应的点值，再由此转化为 \(H\) 的系数表达。
和原来的多项式暴力乘起来相比，两个数的相乘显然轻松多了。
而我们将要学习的 \(DFT\) 即为 将系数表达转化为点值表达，
而 \(IDFT\) 即为 将点值表达转化为系数表达。（\(DFT\) 的逆运算）
DFT 思想
既然我们将 \(F\) 看作一个函数，那么我们将以 \(F(x)\) 来表示往 \(F\) 函数里代入 \(x\) 的点值。
一般来说，\(x\) 显然可以随便取，只要凑够了个数都没问题。
但是我们要考虑 计算的时间复杂度。
如果我们随便取 \(x\) ，那么暴力计算 \(F\) 的一个点值的复杂度就有 \(O(n)\)。
何况我们还要凑够足够多的点值。
如此庞大的复杂度，显然我们承受不起。
所以我们必须从代入的值入手，慎重考虑 \(x\) 的取值。
复数
首先要知道一个 复数 是由它的 实部 和 虚部 组成的。
形如 \(a+bi\) 的表示，其中 \(a\) 即为该数的 实部 ，\(b\) 即为该数的 虚部 。
其中 \(i^2=-1\).
我们可以建立一个二维坐标系（即复平面），\((a,b)\) 即为复数 \(a+bi\) 的坐标。
复平面 中，x轴叫做 实轴 ，y轴叫做 虚轴。
不想画图，脑补一下。
与此同时，复数之间也可以进行运算。
-
复数相加：将对应的 实部 与 虚部 分别相加即可，如 \((a+bi)+(c+di)=(a+c)+(b+d)i\)
-
复数相减：将对应的 实部 与 虚部 分别相减即可，如 \((a+bi)-(c+di)=(a-c)+(b-d)i\)
-
复数相乘：像多项式一样相乘即可，如 \((a+bi) \times (c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i\)
-
复数共轭：实部相同，虚部相反，如 \(a+bi\) 与 \(a-bi\) 共轭。
特殊的，\((a+bi) \times (a-bi)=a^2+b^2\) (为 实数) -
复数相除：分子分母 同乘 分母的共轭 使得分母有理化，之后自行化简即可，如 \(\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i\)
其中，我们重点探讨一下 复数相乘 在 复平面 上的表现。
首先明确两个概念：模长 和 辐角。
-
模长：在 复平面 上连接 原点 和 复数所表示的点 的线段的长度。
-
辐角：在 复平面 上连接 原点 和 复数所表示的点 的 线段 与 实轴 向正方向所形成的夹角。
定理：在两复数相乘时，模长相乘，辐角相加。
证明：
- 模长相乘：设相乘的两个 复数 在 复平面 上对应的点为 \(A\) 和 \(B\).
其中，\(A\) 对应的复数为 \(a+bi\)，\(B\) 的为 \(c+di\).
则它们相乘得到的复数为 \((ac-bd)+(ad+bc)i\) ，设其在复平面上对应的点为 \(C\).
\(A\) 的坐标为 \((a,b)\) ，根据 勾股定理 计算其模长为 \(\sqrt{a^2+b^2}\).
\(B\) 的坐标为 \((c,d)\) ，同理计算其模长为 \(\sqrt{c^2+d^2}\).
\(C\) 的坐标为 \((ac-bd,ad+bc)\) ，计算其模长为
\(\sqrt{(ac-bd)^2+(ad+bc)^2}=\sqrt{a^2c^2-2abcd+b^2d^2+a^2d^2+2abcd+b^2c^2}\)
\(=\sqrt{a^2c^2+a^2d^2+b^2c^2+b^2d^2}\)
\(=\sqrt{(a^2+b^2)(c^2+d^2)}\)
\(=\sqrt{a^2+b^2} \sqrt{c^2+d^2}\)
\(=AO \cdot BO\)
证毕。
- 辐角相加：证相似三角形即可。
太麻烦不想写
复数运算参考代码：
struct CP {
double x,y;
friend CP operator + (const CP& a,const CP& b)
{return (CP){a.x+b.x,a.y+b.y};}
friend CP operator - (const CP& a,const CP& b)
{return (CP){a.x-b.x,a.y-b.y};}
friend CP operator * (const CP& a,const CP& b)
{return (CP){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
};
单位根
单位根是 n次幂为1的复数
换句话说，即 \(x^n=1\) 的复数解。
它的模长也显然等于 \(1\).
我们考虑复平面上的 单位圆 ，即以原点为圆心，半径为 1 的圆。
考虑一个复数，无非有三种情况：
-
若 \(|x|>1\) ，则 \(|x|^n>1\) ，即越乘越大
-
若 \(|x|<1\) ，则 \(|x|^n<1\) ，即越乘越小
-
若 \(|x|=1\) ，则 \(|x|^n=1\)
综上所述，仅当复数 \(x\) 的模长等于 \(1\) 时，才有可能是 单位根 。
而模长为 \(1\) 的点全部都在单位圆上。
进而推得，单位根都在单位圆上。
我们用符号 \(\omega\) 来表示单位根。
而 \(\omega_{n}^{k}\) 则表示 \(n\) 次单位根的 \(k\) 次幂。
显然 \(\omega_n^0=1\) ，而 \(\omega^1_n\) 为其逆时针旋转 \(\frac{1}{n}\) 圆，\(\omega_n^2 … \omega_n^{n-1}\) 以此类推。
至于计算 \(\omega^k_n\) 的坐标，我们可以利用之前学 三角函数 时的知识来求。
根据上述意义，容易得到 \(\omega_n^k (\ cos(\frac{k}{n} \cdot 2\pi)\ ,\ sin(\frac{k}{n} \cdot 2\pi)\ )\) .
在 \(k<0\) 或 \(n \le k\) 时，将 \(k\) 模上一个 \(n\) 即可。
通过这些方法，我们就可以计算出单位根的坐标。
快速DFT
我们尝试以分治形式来实现 DFT .
考虑将 \(F\) 分成左右两半：一半偶次幂，一半奇次幂。
为了满足分治的条件，我们再将其转换成相同的形式。
即：\(F_L=F[0]+F[2]x+F[4]x^2+…+F[n-2]x^{n/2-1}\)
\(F_R=F[1]+F[3]x+F[5]x^2+…+F[n-1]x^{x/2-1}\)
容易知道：\(F(x)=F_L(x^2)+xF_R(x^2)\) (tips：这是点值)
接下来，我们要把刚刚学习的单位根代入到里面，从而产生一些神奇的事情。
我们先令 \(k<n/2\).
-
当 \(m < n/2\) 时，\(F(\omega^m_n)=F(\omega^k_n)=F_L(\omega^k_{n/2})+\omega_n^k F_R(\omega^k_{n/2})\)
-
当 \(m \ge n/2\) 时，\(F(\omega^m_n)=F(\omega^{k+n/2}_n)=F_L(\omega^k_{n/2})-\omega^k_nF_R(\omega^k_{n/2})\).
至于为什么要化成差不多的形式？
不当冤大头，省下求不必要的值的时间。
Tips：
1.\(\omega^{k+n/2}_n=-\omega^k_n\) ，（相当多转了半圈，与原来位置相反）
2.\(\omega^{2k}_n=\omega^k_{n/2}\) ，（相当于原来的每份翻了一倍，即份数少了一半）
可以画画图理解一下。
可以发现，在代入了这些单位根后，我们的转移式也恰好满足了分治的条件。
所以我们只要一直分治下去，就能够求出点值了。
慢着！
分治的转移每次都是对半分，万一 \(n\) 不是 \(2\) 的整幂怎么办？
补 前导 \(0\) 即可，既可以补充位置，也不影响我们的计算。
了解这些，我们就可以开始写 DFT 了。
参考代码：
void DFT(vector<CP>&F) { //养成卡常好习惯
int n=F.size();
if(n==1) return ;
vector<CP>FL,FR;
FL.resize(n/2);
FR.resize(n/2);
for(int i=0;i<n/2;i++)
FL[i]=F[i<<1],FR[i]=F[i<<1|1];
DFT(FL);
DFT(FR);
int k=0;
while((1<<k)<n) k++;
for(int i=0;i<n/2;i++) {
CP tmp=w[k][i]*FR[i];
F[i]=FL[i]+tmp;
F[i+n/2]=FL[i]-tmp;
}
}
在求出 \(F,G\) 的点值后，别忘了将它们对应点值相乘得到 \(H\) 的点值。
参考代码：
for(int i=0;i<len;i++)
F[i]=F[i]*G[i];
IDFT 思想
在学会 DFT 后，我们只得到了 \(H\) 的点值表达，我们还需要应题目要求，将其转化回 系数表达。
既然 DFT 是将 系数表达转化为点值表达 ，那我们还要学会它的逆运算， 将点值表达转化为系数表达，也即 IDFT.
现在我们手上有 \(H\) 的点值表达，不妨设它们为 数列 \(h\).
显然，\(h[k]=H(\omega^k_n)=\sum \limits_{i=0}^{n-1}H[i] \cdot (\omega^k_n)^i\).
而由这个式子，我们可以推得：\(n \cdot H[k]=\sum \limits_{i=0}^{n-1}h[i] \cdot (\omega^{-k}_n)^i\).
证明：
然后我们分类讨论：
-
若 \(j=k\) ，则 贡献 \(=n \cdot H[k]\) ，（有 \(n\) 个符合要求的 \(j\) ）
-
若 \(j \neq k\) ，设 \(p=j-k\) ，则 贡献 \(=\sum \sum \limits_{i=0}^{n-1} \omega^{ip}_n H[k+p]=\sum \omega^p_n \cdot H[k+p] \sum \limits_{i=0}^{n-1} \omega^i_n\)
我们用等比数列求和公式可以得到：
\(\sum \limits_{i=0}^{n-1} \omega^i_n = \frac{\omega^{np}_n-1}{\omega^p_n-1}=\frac{\omega^0_n-1}{\omega^p_n-1}=\frac{1-1}{\omega^p_n-1}=0\)
我们再将其代回贡献，可得这部分贡献为 \(0\).
综上，总和即为 \(n \cdot H[k]\).
证毕。
然后我们就可以利用 \(n \cdot H[k]=\sum \limits_{i=0}^{n-1}h[i] \cdot (\omega^{-k}_n)^i\) 来进行 IDFT 了。
我们还可以发现，\(\omega^{-i}_n = \omega^{n-i}_n\) .
只要将 \(h\) 对应区间翻转一下，即可对应相应位置，从而又省下一堆时间。
除此之外，可以发现 IDFT 中部分运算与 DFT 相似，合理利用即可。
废话不多说，看代码：
void IDFT(vector<CP>&F) {
int n=F.size();
DFT(F);
reverse(&F[1],&F[n]);
for(int i=0;i<n;i++)
F[i].x/=n,F[i].y/=n;
}
至此，我们已经完成了 \(FFT\) 的重要步骤学习了。
时间复杂度：\(O(nlog(n))\).
当然，要记得卡常。
最后附上刚好卡过模板题的代码：
//FFT
#include<bits/stdc++.h>
using namespace std;
const int maxn=1e6+5;
const double Pi=acos(-1);
struct CP {
double x,y;
friend CP operator + (const CP& a,const CP& b)
{return (CP){a.x+b.x,a.y+b.y};}
friend CP operator - (const CP& a,const CP& b)
{return (CP){a.x-b.x,a.y-b.y};}
friend CP operator * (const CP& a,const CP& b)
{return (CP){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
};
vector<CP>w[20];
void Init(int N) {
for(int k=1;(1<<k)<=N;k++) {
int n=1<<k;
w[k].resize(n/2);
for(int i=0;i<n/2;i++)
w[k][i]=(CP){cos(i*2*Pi/n),sin(i*2*Pi/n)};
}
}
void DFT(vector<CP>&F) {
int n=F.size();
if(n==1) return ;
vector<CP>FL,FR;
FL.resize(n/2);
FR.resize(n/2);
for(int i=0;i<n/2;i++)
FL[i]=F[i<<1],FR[i]=F[i<<1|1];
DFT(FL);
DFT(FR);
int k=0;
while((1<<k)<n) k++;
for(int i=0;i<n/2;i++) {
CP tmp=w[k][i]*FR[i];
F[i]=FL[i]+tmp;
F[i+n/2]=FL[i]-tmp;
}
}
void IDFT(vector<CP>&F) {
int n=F.size();
DFT(F);
reverse(&F[1],&F[n]);
for(int i=0;i<n;i++)
F[i].x/=n,F[i].y/=n;
}
int main() {
int n,m,len=1;
scanf("%d%d",&n,&m);
while(len<=n+m) len<<=1;
Init(len);
vector<CP> F,G;
F.resize(len),G.resize(len);
for(int i=0;i<=n;i++)
scanf("%lf",&F[i].x);
for(int i=0;i<=m;i++)
scanf("%lf",&G[i].x);
DFT(F);
DFT(G);
for(int i=0;i<len;i++)
F[i]=F[i]*G[i];
IDFT(F);
for(int i=0;i<=n+m;i++) {
int tmp=F[i].x+0.1;
printf("%d ",tmp);
}
return 0;
}
浙公网安备 33010602011771号