FFT学习笔记
$\quad $ 本人蒟蒻，只能介绍FFT在OI中的应用，如有错误或不当之处还请指出。
$\quad $ 首先先说一下那一堆什么什么 \(TT\) 的都是什么
DFT：离散傅里叶变换
FFT ：快速傅里叶变换 用于求多项式乘法 \(O(n log(n))\)
FNTT/NTT ：FTT的优化，常数及精度更优
FWT：快速沃尔什变换
MTT：任意模数NTT
FMT 快速莫比乌斯变换
SFT 稀疏傅里叶变换 百度百科说比FFT快10到100倍但是在OI里真的有应用吗？
$\quad $ 后面三个我都不会，这里也不介绍 NTT我也不会，仙姑了~
以下是一些前置知识（当然不只是前置知识🌚）
多项式
$\quad $ 多项式有系数表示法和点值表示法两种，我们常见的 \(1+x+x^2+x^3\) 这样的就叫做系数表示法，而用若干个点表示多项式的就叫做点值表示法。
$\quad $ 把多项式看作函数，如果你学过一点插值的话，就可以知道给定坐标系中横坐标互不相同的 \(n+1\) 个点，可以唯一确定一条 \(n\) 次函数经过这 \(n+1\) 个点。
$\quad $ 用系数表示法来计算两个多项式之积，枚举各项系数相乘，时间复杂度是 \(O(n^2)\) 的。但是点值表示法呢？我们只需要把在某一点处两个多项式的值相乘，就能得到结果多项式在该处的值，直接枚举 \(n\) 个点即可，时间复杂度 \(O(n)\) 。
$\quad $ 那么我们就可以考虑将两个多项式全部转化为点值表示，然后计算乘积，再将结果多项式从点值表示转化为系数表示，这也就是FFT的大致流程。
$\quad $ 那么点值表示和系数表示相互转化的时间复杂度呢？发现是 \(O(n^2)\) 的，《这么一整直接把常数优化大了🌚》
$\quad $ 那么我们就要引出快速傅里叶变换了，在这之前我们引出单位根这个概念。
单位根
$\quad $ 如果有没学过复数运算的，请先去学复数运算。
$\quad $ 首先对于方程 \(x^n=1\) ，其在复数域中有 \(n\) 个解，并且如果画出来的话，这 \(n\) 个解恰好等分单位圆，我们按照逆时针顺序分别叫他们 \(\omega _{n}^{1} 、 \omega _{n}^{2} 、\cdots 、\omega _{n}^{n}\) ，不难发现其性质：
$\quad $ 前三条性质显然，我证明一下第四条：
$\quad $ 首先看这个级数：\(\sum _{i=0}^{n}x^i\) ，我们令 \(f=\sum _{i=0}^{n}x^i\) ，那么 \(xf=\sum _{i=1}^{n}x^i\) ，相减可以求出 \(f=\frac{x ^{n+1}-1}{x-1}\) ，当然在分母为 \(0\) 时不能用该式求解。我们发现上式就是一个等比数列的形式，所以：
\begin{aligned}
\sum _{i=0}^{n-1}{ ( \omega _{n}^{k} ) ^ i}&=\frac{(\omega _{n}^{k}) ^ n}{\omega _{n}^{k}-1}=\frac{(\omega _{n}^{n}) ^ k -1}{\omega _{n}^{k}-1}=\frac{1-1}{\omega _{n}^{k}-1}
\end{aligned}
$\quad $ 这里分子一定为零了，我们讨论分母。发现当且仅当 \(k=n\) 时，分母为零，此时直接代入原式得到答案为 \(n\) 。其余情况该式均为零。
快速傅里叶变换
$\quad $ 那么介绍这个单位根干什么呢？我们可以取一个 \(n-1\) 次多项式在这 \(n\) 个单位根处的取值来作为这个多项式的点值表示，至于为什么这么取，看下面就知道了。
$\quad $ 我们定义一个 \(n-1\) 次多项式 \(A(x)=\sum _{i=0}^{n-1}a _{i}x ^{i}\) ，其各项系数为 \(a _{0},a _{1},\cdots,a _{n-1}\) 。
$\quad $ 我们按照各项系数下标的奇偶性把系数分成两类，并各作为一个 \(\frac{n}{2}-1\) 次多项式的系数，从而得到下面两个新的多项式：
$\quad $ 并且有：\(A(x)=A _{1}(x^2)+xA _{2}(x^2)\)
$\quad $ 再看一下我们要求什么，我们要求：\(A(\omega _{n}^{1}),A(\omega _{n}^{2}),\cdots,A(\omega _{n}^{n})\) 。
$\quad $ 分别代入 \(A(\omega _{n}^{k})\) 和 \(A(\omega _{n}^{k+\frac{n}{2}})\) ，我们可以得到以下两个式子：
\begin{aligned}
A(\omega _{n}^{k})&=A _{1} ((\omega _{n}^{k}) ^ 2)+\omega _{n}^{k} A _{2} ((\omega _{n}^{k}) ^ 2)\\
&=A _{1} (\omega _{n}^{2k})+\omega _{n}^{k} A _{2} (\omega _{n}^{2k})\\
&=A _{1} (\omega _{\frac{n}{2}}^{k})+\omega _{n}^{k} A _{2} (\omega _{\frac{n}{2}}^{k})\\
\end{aligned}
\begin{aligned}
A(\omega _{n}^{k+\frac{n}{2}})&=A _{1} ((\omega _{n}^{k+\frac{n}{2}}) ^ 2)+\omega _{n}^{k+\frac{n}{2}} A _{2} ((\omega _{n}^{k+\frac{n}{2}}) ^ 2)\\
&=A _{1} (\omega _{n}^{2k+n})+\omega _{n}^{k+\frac{n}{2}} A _{2} (\omega _{n}^{2k+n})\\
&=A _{1} (\omega _{\frac{n}{2}}^{k})-\omega _{n}^{k} A _{2} (\omega _{\frac{n}{2}}^{k})\\
\end{aligned}
$\quad $ 发现这两个式子惊人的相似，那我们只需要求出 \(A(\omega _{n}^{1}),A(\omega _{n}^{2}),\cdots,A(\omega _{n}^{\frac{n}{2}})\) ，然后就可以 \(O(1)\) 地求出 \(A(\omega _{n}^{\frac{n}{2}+1}),A(\omega _{n}^{\frac{n}{2}+2}),\cdots,A(\omega _{n}^{n})\) 。
$\quad $ 我们考虑等式右边的东西，它的形式和原来那个问题一样，于是我们就可以递归求解了。但是这里的除法均不是向下取整，所以在实现的时候我们需要把 \(n\) 补成一个 \(2\) 的整数次幂，也就是把原先的多项式都补成一个 \(2^k\) 次的多项式，至于多出来的系数，直接置为零就好了，也不会影响答案。不难发现，这样的算法时间复杂度是 \(O(n log(n))\) 的。
放个码
struct complex{
double x,y;
complex operator+(const complex op)const{
return complex{x+op.x,y+op.y};
}
complex operator-(const complex op)const{
return complex{x-op.x,y-op.y};
}
complex operator*(const complex op)const{
return complex{x*op.x-y*op.y,x*op.y+y*op.x};
}
};//当然你也可以直接用std里的complex
void FFT(complex*a,int n,int op){
if(!n)return;
complex a1[n],a2[n];
for(int i=0;i<n;++i)a1[i]=a[i<<1],a2[i]=a[i<<1|1];
FFT(a1,n>>1,op);FFT(a2,n>>1,op);
complex W{cos(pi/n),sin(pi/n)*op},w{1,0};
for(int i=0;i<n;++i,w=w*W)a[i]=a1[i]+w*a2[i],a[i+n]=a1[i]-w*a2[i];
}
$\quad $ 于是我们有了 \(O(nlog(n))\) 将一个系数表示的多项式转为点值表示的方法，然后呢？不要忘了我们还要将点值表示的结果多项式再转为系数表示，所以下面我们介绍快速傅里叶逆变换。
快速傅里叶逆变换
$\quad $ 这里有线性代数的推法，如果你是会范德蒙德矩阵巨佬的话，就不要看这里的介绍了。🌚
$\quad $ 我们再定义一个多项式 \(B=\sum _{i=0}^{n-1}y _{i} x^i\) ，这里的 \(y _{i}\) 是 \(A\) 在 \(\omega _{n}^{i}\) 处的值，即 \(A(\omega _{n}^{i})\) （说了句废话🌚）。
$\quad $ 上面的 \(\omega _{n}^{0}\) 就等于 \(\omega _{n}^{n}\) ，怕有人没想过来还是写一下🌚。
$\quad $ 这是系数表示 \(B\) ，我们考虑用点值表示法来表示 \(B\) ，这次我们取 \(b _{i}=\omega _{n}^{-i}\) ，那为什么要这么取呢？我也不知道 可以去学一下用范德蒙德矩阵推导，然后就知道为什么是这个了🌚。
$\quad $ 我们先试着把 \(b _{k}\) 代入 \(B\) 中：
\begin{aligned}
B(b _{k})&=\sum _{i=0}^{n-1}A(\omega _{n}^{i})(\omega _{n}^{-k}) ^{i}\\
&=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}a _{j}(\omega _{n}^{i}) ^{j}(\omega _{n}^{-k}) ^{i}\\
&=\sum _{j=0}^{n-1}a _{j}\sum _{i=0}^{n-1}(\omega _{n}^{j-k}) ^{i}
\end{aligned}
$\quad $ 后面的部分是个等比数列，直接用上文中单位根的第四条性质可以得到只有当 \(j=k\) 时，该部分值为 \(n\) ，否则为零，不难得出：
$\quad $ 发现我们只需要得到 \(B\) 在 \(\omega _{n}^{-1},\omega _{n}^{-2},\cdots,\omega _{n}^{-n},\) 处的取值，就可以推出结果多项式 \(A\) 的各项系数，也就有了我们想要的东西。
$\quad $ 如果你有过人的观察能力，就发现这步和上一步相比多了个负号，那么我们是否还可以用之前的方法去做呢？当然是可以的，这里我就不再赘述了，自己推导即可。实现的时候多传一个参数区分开变换和逆变换就行了。
$\quad $ 于是我们就得到了递归版的FFT：
点击查看代码
#include<bits/stdc++.h>
const int N=4e6+10;
const double pi=acos(-1);
struct complex{
double x,y;
complex operator+(const complex op)const{
return complex{x+op.x,y+op.y};
}
complex operator-(const complex op)const{
return complex{x-op.x,y-op.y};
}
complex operator*(const complex op)const{
return complex{x*op.x-y*op.y,x*op.y+y*op.x};
}
}a[N],b[N];
int n,m;
inline int read(){
int ans=0;char ch=getchar_unlocked();bool fl=0;
while(ch<'0'||ch>'9'){if(ch=='-')fl=1;ch=getchar_unlocked();}
while(ch>='0'&&ch<='9')ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar_unlocked();
return fl?~ans+1:ans;
}
inline void print(int x){(x<0)&&(putchar_unlocked('-'),x=-x);(x>9)&&(print(x/10),0);putchar_unlocked(x%10|48);}
void FFT(complex*a,int n,int op){
if(!n)return;
complex a0[n],a1[n];
for(int i=0;i<n;++i)a0[i]=a[i<<1],a1[i]=a[i<<1|1];
FFT(a0,n>>1,op);FFT(a1,n>>1,op);
complex W{cos(pi/n),sin(pi/n)*op},w{1,0};
for(int i=0;i<n;++i,w=w*W)a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}
int main(){
// freopen("1.in","r",stdin);
scanf("%d%d",&n,&m);
for(int i=0;i<=n;++i)a[i]={1.0*read(),0};
for(int i=0;i<=m;++i)b[i]={1.0*read(),0};
m=m+n;n=1;
while(n<=m)n<<=1;
FFT(a,n>>1,1);FFT(b,n>>1,1);
for(int i=0;i<n;++i)a[i]=a[i]*b[i];
FFT(a,n>>1,-1);
for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
}
$\quad $ 这里应该是过不了板子题的，但是它意外过了，没看懂，但是还是要引出优化。《就当上面那个码过不了吧》
$\quad $ 于是我们引出一些优化，注意到上面的代码中有：a[i]=a1[i]+w*a2[i],a[i+n]=a1[i]-w*a2[i]; ，其中 w*a2[i] 我们计算了两次，而复数相乘常数会很大，所以我们拿个变量记一下就好了。于是我们完成了所有的优化
$\quad $ 《显然不是的🌚》，发现递归常数大，我们考虑是否可以摆脱递归，就有了下面这个优化：
位逆序置换
$\quad $ 我们观察每次分治后下标的位置变化，这里以七次多项式举例：
\([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}]\)
\([a _{0}][a _{4}][a _{2}][a _{6}][a _{1}][a _{5}][a _{3}][a _{7}]\)
$\quad $ 再用二进制表示出开始和结尾的位置：
\(000,001,010,011,100,101,110,111\)
\(000,100,010,110,001,101,011,111\)
$\quad $ 如果我们知道每个系数递归到最后的位置，那么我们就可以自下而上的递推求解，也就摆脱了递归。
$\quad $ 那么我们如何求呢？显然我们可以一位一位地翻转求解，这样的时间复杂度是 \(O(n log(n))\) 的，当然我们还有更好的方法。
$\quad $ 记原先下标为 \(n\) 的系数最后的位置为 \(f(n)\) ，我们考虑递推求解：
$\quad $ 我们可以先将 \(n\) 右移一位，翻过来再右移一位，最后再加上 \(n\) 最低位翻转的结果，即可得到 \(f(n)\) ，写成数学式子是：\(f(n)=\lfloor \frac{f(\lfloor\frac{n}{2}\rfloor)}{2}\rfloor+(n\bmod2)\times \frac{len}{2}\) 。证明我也不知道
非递归版FFT
#include<bits/stdc++.h>
const int N=4e6+10;
const double pi=acos(-1);
struct complex{
double x,y;
complex operator+(const complex op)const{
return complex{x+op.x,y+op.y};
}
complex operator-(const complex op)const{
return complex{x-op.x,y-op.y};
}
complex operator*(const complex op)const{
return complex{x*op.x-y*op.y,x*op.y+y*op.x};
}
}a[N],b[N],buf[N];
int n,m,pos[N];
inline int read(){
int ans=0;char ch=getchar_unlocked();bool fl=0;
while(ch<'0'||ch>'9'){if(ch=='-')fl=1;ch=getchar_unlocked();}
while(ch>='0'&&ch<='9')ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar_unlocked();
return fl?~ans+1:ans;
}
inline void print(int x){(x<0)&&(putchar_unlocked('-'),x=-x);(x>9)&&(print(x/10),0);putchar_unlocked(x%10|48);}
void FFT(complex a[],int flag){
for(int i=0;i<n;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]);
for(int mid=1;mid<n;mid<<=1){
complex W={cos(pi/mid),flag*sin(pi/mid)};
for(int r=mid<<1,j=0;j<n;j+=r){
complex w={1,0};
for(int k=0;k<mid;k++,w=w*W){
complex op=w*a[j+k+mid];
buf[j+k]=a[j+k]+op;
buf[j+k+mid]=a[j+k]-op;
}
}
for(int i=0;i<n;i++)a[i]=buf[i];
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<=n;++i)a[i]={1.0*read(),0};
for(int i=0;i<=m;++i)b[i]={1.0*read(),0};
m=m+n;n=1;
while(n<=m)n<<=1;
for(int i=1;i<=n;i++)pos[i]=(pos[i>>1]>>1)|((i&1)*(n>>1));
FFT(a,1);FFT(b,1);
for(int i=0;i<n;++i)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
}
$\quad $ 发现递归改成非递归，直接慢了0.49s，这是怎么烩逝呢？我也不知道，但理论上是快了🌚，大概是我写的常数太大了吧🌚。
$\quad $ 但是无所谓，我们认为它起到了一定的优化作用，但是我们发现多了一个数组，那我们可不可以不要这个数组呢？于是就有了下面这个优化：
蝶形运算优化
$\quad $ 原本的方法可能要动脑子，我这里就说一种不用动脑子的方法。
$\quad $ 观察这里：
for(int r=mid<<1,j=0;j<n;j+=r){
complex w={1,0};
for(int k=0;k<mid;k++,w=w*W){
complex op=w*a[j+k+mid];
buf[j+k]=a[j+k]+op;
buf[j+k+mid]=a[j+k]-op;
}
}
for(int i=0;i<n;i++)a[i]=buf[i];
$\quad $ 发现两次赋值比较慢，直接拿两个变量临时存一下，直赋值给原数组即可，这样就少了个数组，并且不会有影响。
就是这样
#include<bits/stdc++.h>
const int N=4e6+10;
const double pi=acos(-1);
struct complex{
double x,y;
complex operator+(const complex op)const{
return complex{x+op.x,y+op.y};
}
complex operator-(const complex op)const{
return complex{x-op.x,y-op.y};
}
complex operator*(const complex op)const{
return complex{x*op.x-y*op.y,x*op.y+y*op.x};
}
}a[N],b[N];
int n,m,pos[N];
inline int read(){
int ans=0;char ch=getchar_unlocked();bool fl=0;
while(ch<'0'||ch>'9'){if(ch=='-')fl=1;ch=getchar_unlocked();}
while(ch>='0'&&ch<='9')ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar_unlocked();
return fl?~ans+1:ans;
}
inline void print(int x){(x<0)&&(putchar_unlocked('-'),x=-x);(x>9)&&(print(x/10),0);putchar_unlocked(x%10|48);}
void FFT(complex a[],int flag){
for(int i=0;i<n;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]);
for(int mid=1;mid<n;mid<<=1){
complex W={cos(pi/mid),flag*sin(pi/mid)};
for(int r=mid<<1,j=0;j<n;j+=r){
complex w={1,0};
for(int k=0;k<mid;k++,w=w*W){
complex op=w*a[j+k+mid],opl=a[j+k];
a[j+k]=opl+op;
a[j+k+mid]=opl-op;
}
}
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<=n;++i)a[i]={1.0*read(),0};
for(int i=0;i<=m;++i)b[i]={1.0*read(),0};
m=m+n;n=1;
while(n<=m)n<<=1;
for(int i=1;i<=n;i++)pos[i]=(pos[i>>1]>>1)|((i&1)*(n>>1));
FFT(a,1);FFT(b,1);
for(int i=0;i<n;++i)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
}
$\quad $ 一个挺简单也挺好的优化，比着上一版非递归直接少了1.07s，想再了解点的移步oi-wiki，当然还有一个优化
三步变两步优化
$\quad $ 顾名思义，之前的代码调用了三次FFT函数，我们实际上是可以只调用两次的。
$\quad $ 对于给定的两个多项式 \(A(x)=\sum _{i=0}^{n-1}a _{i}x ^{i},B(x)=\sum _{i=0}^{m-1}b _{i} x^i\) ，我们令 \(a、b\) 分别为新的多项式各项的虚部和实部（我这里用 \(\rm i\) 来表示虚数单位，注意辨别），从而得到多项式 \(H(x)=\sum _{i=0}^{len-1}(a _{i} + b _{i} {\rm i}) x ^i\)
$\quad $ 可以将其拆成两项相乘：
$\quad $ 这样我们只需要计算 \(H^2\) ，再取其虚部的 \(\frac{1}{2}\) 即可，但是这样做的误差会大一点。
$\quad $ 到这里FFT的内容就结束了（其实还有一个DIF DIT来着，但我还没学会）。
后记
$\quad $ 开始看FFT的时候感觉好难，甚至以为单位根是最难的东西，后来发现这比着后面的部分简单的要多。许多东西也是边学边写的，如果有不当之处还请各位大佬指出，本人将不胜感激。
$\quad $ 最后再写个无关的东西放松大脑吧
三项式定理
$\quad $ 其实这个是可以推广为 \(n\) 项式的。自己探索吧（其实是我觉得式子太麻烦不想写了）。
