快速傅里叶变换（FFT）学习笔记
前言
本文仅作为本人学习笔记使用。但若您的实力比我强，当然可以阅读。
基础知识
FFT，Fast Fourier Transform，快速傅里叶变换。
是一种值域转频域的牛逼变换 （没用），在 OI 中常用的形式是 DFT+IDFT，用于 \(O(n\log n)\) 计算多项式乘法（或卷积）。
记号
- \(F(x)\)：多项式 \(F\)，如 \(F(x)=x^2-x+2\)。
- 项数，\(n\)。对于 \(F(x)=x^2-x+2\)，\(n=3\)，叫它 \(n-1=2\) 次多项式。
- 系数。系数都是参数（非未知数），记 \(F[i]\) 为 \(F(x)\) 的 \(i\) 次项系数。
对于 \(F(x)=x^2-x+2\)，\(F[2]=1,F[1]=-1,F[0]=2\)。 - 卷积，记做 \(*\)。即\[\large H(x)=F(x)*G(x) \]对于乘法，它的卷积形式（仅系数）为\[\large H[k]=\sum_{\mathclap{i+j=k}}F[i]\times G[j] \]
点值表达
对于一个二次多项式 \(F(x)=ax+b\)，它表达了平面直角坐标系上的一条直线，如果给定两个这条线上的点，显然可以两点确定一直线，求出解析式，即 \(a,b\) 的值。
更一般地，对于一个 \(n\) 次多项式，只需 \(n+1\) 个点，就可以确定它们所在的多项式。
于是我们就可以将多项式 \(F(x)\) 表达为 \(n+1\) 个点，我们叫它点值表达法。
记 \(F(x)\) 的点值表达法为 \(\large (FX,FY)=(\{Fx_0,\cdots,Fx_n\},\{Fy_0,\cdots,Fy_n\})\)，满足 \(F(Fx_i)=Fy_i\)。
对 \(F(x)\) 和 \(G(x)\) 的暴力卷积显然是 \(O(n^2)\) 的，但转化为点值表达法后只需 \(O(n)\)，前提是点的横坐标一一对应。
具体地，如果我们知道了对于 \(n\) 个点的横坐标集 \(X=\{x_0,x_1,\cdots,x_{n-1},x_n\}\) 中每个 \(x_i\) 对应的纵坐标值 \(F(x_i)=Fy_i,G(x_i)=Gy_i\)，那么它们的乘积 \(H(x)\) 的点值表达中显然必然有 \(HY=\{Fy_1Gy_1,\cdots,Fy_nGy_n\}\)。
但相乘后，\(H(x)=F(x)G(x)\) 为 \(2n\) 次多项式，需要 \(2n+1\) 个点确定，所以需要选取 \(HX,|HX|=2n+1\) 个横坐标计算 \(HY\)。
点值表达法好处：
- 相乘快。
- 自选 \(X\)。
缺漏：
- 相乘后为点值表达法，难以得知解析式。
- 不知道如何选择 \(HX\)。
FFT 思想：取 \(HX\)，计算 \(FY,GY\) \(\to\) 计算 \(HY\) \(\to\) 得到解析式（难）。
简称：优点就是缺点。所以我们要解决问题，当然得从选择 \(HX\) 入手。一会儿会讲。
复数
DFT有一个妙处,就是代入什么由你自己决定,只要点值个数足够。
我们这些蒟蒻第一感觉都会选择人畜无害的正整数,或者有理数什么的。
但是这些东西虽然在普通人看来计算简单,但并没有什么能够加以利用的优秀性质。
事实证明,找一些毒瘤东西代入进去是个好想法。
——command_block
复数，形如 \(a+b\text{i}\) 的数。在复平面上。\(\text{i}^2=-1\)。
运算：
- 加减：\((a+b\text{i})\pm(c+d\text{i})=(a\pm c)+(b\pm d)\text{i}\)
- 乘：类似一次多项式，\((a+b\text{i})\times(c+d\text{i})=(ac-bd)+(ad+bc)\text{i}\)
- 共轭：\(a+b\text{i}\) 的共轭为 \(a-b\text{i}\)，二者相乘为实数 \(a^2+b^2\)。
- 除：化简时同乘共轭。（有点像根式化简呢）。
- 模：\(a+b\text{i}\) 的模记为 \(|a+b\text{i}|=\sqrt{a^2+b^2}\)。
- 幅角：\(\arctan \cfrac{b}{a}\)。即 始边为 \(x\) 轴，终边为 \(O\) 与 \(a+b\text{i}\) 所表示的点的间的直线 的夹角。
结论：两复数相乘，模长相乘，幅角相加。
特别地，如果模长恒为 \(1\)（村规：叫它单位圆复数，顾名思义复平面上的单位圆上的复数），即特化为：两单位圆复数相乘，幅角相加。（做铺垫。）
单位根
简单说，\(n\) 次单位根就是方程
的复数解（定义）。
显然，只有单位圆复数才能成为单位根。
理由如下：若复数 \(x\)，
- 若 \(|x|>1,|x^n|=|x|^n>1\)，它的 \(n\) 次幂会向外远离单位圆。
- 若 \(|x|<1,|x^n|=|x|^n<1\)，它的 \(n\) 次幂会向内远离单位圆。
- 若 \(|x|=1\)，即 \(x\) 为单位圆复数，\(|x^n|=|x|^n=1^n=1\)，它的 \(n\) 次幂永远停留在单位圆上。
易得：幅角为 \(\large 0,\cfrac{1}{n}\times 2\pi,\cdots,\cfrac{n-1}{n}\times 2\pi\) 的单位圆复数为 \(n\) 次单位根。
即 \(n\) 次单位根有 \(n\) 个，它们 \(n\) 等分圆。
幅角为 \(0,\cfrac{1}{n}\times 2\pi,\cdots,\cfrac{n-1}{n}\times 2\pi\) 的 \(n\) 次单位根依次记为：\(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\)（记号）。
\(n\) 次单位根的上标其实就是幂，下面的 性质 可以助你理解：
- \(\omega_n^0=1\)：显然的，为方程 \((1)\) 的实数解。同时也可以理解为除 \(0\) 外任何数的 \(0\) 次方都为 \(1\)。
- \(\omega_n^k=\left(\omega_n^1\right)^k\)：\(\omega_n^1\) 幅角为 \(\cfrac{1}{n}\times 2\pi\)，自乘 \(k\) 次后，根据“两单位圆复数相乘，幅角相加”，幅角为 \(\cfrac{k}{n}\times 2\pi\)，所以为 \(\omega_n^k\)。当然直接用幂理解也是可以的。任何一个 \(n\) 次单位根都可以写成 \(\omega_n^1\) 的整幂。
- \(\omega_n^j\times \omega_n^k=\omega_n^{j+k},(\omega_n^j)^k=\omega_n^{jk}\)：类似性质 2，可用幂解释。
- \(\omega_n^{k+n/2}=-\omega_n^k\)：\(\omega_n^{k+n/2}=\omega_n^k\times\omega_n^{n/2}\)，\(\omega_n^{n/2}\) 的幅角为 \(\cfrac{n/2}{n}\times 2\pi=\pi\)，根据“两单位圆复数相乘，幅角相加”，也就是将 \(\omega_n^k\) 旋转 \(180\degree\)。复平面上，为 \(-\omega_n^k\)。
- \(\left(\omega_n^k\right)^2=\omega_{n/2}^k\)：\(2\times\cfrac{1}{n}=\cfrac{2}{n}\)，自己理解吧。
“任何一个 \(n\) 次单位根都可以写成 \(\omega_n^1\) 的整幂”，既然如此，如何算出 \(\omega_n^1\)？
\(\omega_n^1\) 的幅角为 \(\cfrac{1}{n}\times 2\pi=\cfrac{2\pi}{n}\)。根据单位圆上的三角函数
DFT
看清楚，别眨眼。
设 \(n-1\) 次多项式
奇偶分开
保证 \(n\) 为 \(2\) 的整幂，即不会分不均匀。
设两个 \(\cfrac{n}{2}-1\) 次多项式
可得
设 \(0\le k<\cfrac{n}{2}\)，将 \(\omega_n^k\) 代入 \((3)\)，
再将 \(\omega_n^{n/2+k}\) 代入 \((3)\)，
对比 \((4),(5)\) 两式，
肉眼可见，等号右边只有 \(+,-\) 的不同。你这说的屁用没有。
意义：如果我们知道了 \(F_l(x),F_r(x)\) 在 \(\omega_{n/2}^0,\omega_{n/2}^1,\cdots,\omega_{n/2}^{n/2-1}\) 处的点值表示，
套用式子 \((4),(5)\)，就可以分别 \(O(n)\) 求出 \(F(x)\) 在 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n/2-1}\) 处和 \(\omega_n^{n/2},\omega_n^{n/2+1},\cdots,\omega_n^{n-1}\) 处的点值表示。
也就是说，如果我们获得了 \(F(x)\) 的点值表示！
但是有一个啸啸的问题：“如果知道了 \(F_l(x),F_r(x)\) 在 \(\omega_{n/2}^0,\omega_{n/2}^1,\cdots,\omega_{n/2}^{n/2-1}\) 处的点值表示”。我们不知道！
怎么办？暴力代入计算吗？肯定不是。
容易发现，将 \(F(x)\) 的点值表示拆分，通过计算 \(F_l(x),F_r(x)\) 的点值表示来计算，是一个分治的过程。
到最后只剩一项时，\(F(x)\) 为常数，其点值表达式为其本身，直接回溯即可。
小练习：自己推导一遍式 \((4),(5)\)。
这就完了吗？并没有。“保证 \(n\) 为 \(2\) 的整幂，即不会分不均匀。”
没事，大不了将 \(n\) 补位 \(2\) 的整幂，即加几个系数为 \(0\) 的高次项即可。
IDFT
IDFT，实际上就是 DFT 的逆变换。理论也十分简单。
假如将 \(F(x)\) DFT 后得到的纵坐标集为 \(Y\)，那么代入过程为
这里给出 IDFT 的结论：
注意 \(i\)（变量）非 \(\text{i}\)（虚数单位）。
证明：
将 \((6)\) 代入式 \((7)\) 右边
设 \(l=k-i\)，若 \(l=0,k=i\)
若 \(l\ne 0,k=i+l\)
两式相加（合并 \(l=0,l\ne 0\) 两种情况的贡献）
回来看式 \((7)\)
再看看式 \((6)\)
是不是有 十分甚至九分的相似？
我们只需将式 \((6)\) 中 \(\left(\omega_n^j\right)^k\to\left(\omega_n^{-i}\right)^j\)，再 \(/n\) 即可。
也就是说，我们只需把 \(Y\) 当做系数，扔给 DFT（DFT：“我太难了”），再让 DFT 计算时用 \(\omega_n^0,\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-(n-1)}\) 代入即可。它们的计算同 \(\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}\)，用 \(\omega_n^{-k}=\left(\omega_n^{-1}\right)^k\) 计算。
其中
理论存在——总结
实践开始！——代码
递归版
预处理了 \(\omega\)。
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
constexpr int N=1<<21|1;
constexpr double pi=3.141592653589793;
struct cp
{
double a,b;
cp():a(0),b(0){}
cp(double _a,double _b):a(_a),b(_b){}
inline cp operator + (cp y){return cp(a+y.a,b+y.b);}
inline cp operator - (cp y){return cp(a-y.a,b-y.b);}
inline cp operator * (cp y){return cp(a*y.a-b*y.b,a*y.b+b*y.a);}
// inline cp operator ! (){return cp(a,-b);}// gong'e
inline void print(){printf("%lf+%lfi",a,b);}
};
int n,m;
cp w[N];
void init()// calc omega
{
w[n]=w[0]={1,0};
w[1]={cos(2*pi/n),sin(2*pi/n)};
for(int i=2;i<n;i++)w[i]=w[i-1]*w[1];
}
int cn,dis;
cp getw(int i,int base)// get omega
{
dis=cn-__builtin_ctz(base);// n/base=1<<dis
return w[i<<dis];// w_base^i=w_{base<<dis}^{i<<dis}=w_n^{i<<dis}
}// good!
cp f[N],g[N];
int rev[N];
// flg=1:dft;flg=0:idft
void fft(cp *f,int len=n)// dft&idft
{
if(len==1)return;
fft(f,len>>1),fft(f+(len>>1),len>>1);
cp t;
for(int k=0;k<(len>>1);k++)
{
t=getw(k,len)*f[(len>>1)+k];
// printf("w_%d^%d=",len,k),getw(k,len).print(),puts("");
f[(len>>1)+k]=f[k]-t;
f[k]=f[k]+t;
}
}
inline void swp(cp* f){for(int i=0;i<n;i++)if(rev[i]<i)swap(f[rev[i]],f[i]);}
int main()
{
scanf("%d %d",&n,&m);
for(int i=0;i<=n;i++)scanf("%lf",&f[i].a);
for(int i=0;i<=m;i++)scanf("%lf",&g[i].a);
for(m+=n,n=1;n<=m;n<<=1,cn++);// up to 2^cn
for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
swp(f),swp(g),init(),fft(f),fft(g);/* puts(""); */
for(int i=0;i<n;i++)f[i]=f[i]*g[i];
swp(f),reverse(w,w+n+1),fft(f);
for(int i=0;i<=m;i++)printf("%d ",int((f[i].a/n)+.5));
return 0;
}
递推（迭代）版
未预处理 \(\omega\)，更快。
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
constexpr int N=1<<21|1;
constexpr double pi=3.141592653589793;
struct cp
{
double a,b;
cp():a(0),b(0){}
cp(double _a,double _b):a(_a),b(_b){}
inline cp operator + (cp y){return cp(a+y.a,b+y.b);}
inline cp operator - (cp y){return cp(a-y.a,b-y.b);}
inline cp operator * (cp y){return cp(a*y.a-b*y.b,a*y.b+b*y.a);}
// inline void print(){printf("(%lf+%lfi)",a,b);}
};
int n,m;
int rev[N];
void init()// calc omega
{
for(m+=n,n=1;n<=m;n<<=1);
for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
}
void fft(cp *f,int sgn=1)
{
for(int i=0;i<n;i++)if(rev[i]<i)swap(f[rev[i]],f[i]);
cp t,w1,w;
for(int l=2,hl;l<=n;l<<=1)
{
hl=l>>1;
w1={cos(2*pi/l),sgn*sin(2*pi/l)};
for(int i=0;i<n;i+=l)
{
w={1,0};
for(int k=i;k<i+hl;k++)
{
t=w*f[hl+k];
f[hl+k]=f[k]-t;
f[k]=f[k]+t;
w=w*w1;
}
}
}
}
cp f[N],g[N];
int main()
{
scanf("%d %d",&n,&m);
for(int i=0;i<=n;i++)scanf("%lf",&f[i].a);
for(int i=0;i<=m;i++)scanf("%lf",&g[i].a);
init(),fft(f),fft(g);
for(int i=0;i<n;i++)f[i]=f[i]*g[i];
fft(f,-1);
for(int i=0;i<=m;i++)printf("%d ",int(f[i].a/n+.5));
return 0;
}
完结撒花！
全文完。
也许什么时候学一下 NTT？其实不急。
