FFT快速傅里叶变换
DFT：离散傅里叶变换—>\(O(n^2)\)计算多项式乘法
FFT：快速傅里叶变换—>\(O(n\log n)\)计算多项式乘法
FNTT/NTT：快速傅里叶变换的优化版—>优化常数及误差
FWT：快速沃尔什变换—>利用类似FFT的东西解决一类卷积问题
MTT：毛爷爷的FFT—>非常nb/任意模数
FMT 快速莫比乌斯变化
快速傅里叶变换
\(FFT\) 算法的基本思想是分治。就 \(DFT\) 来说，它分治地来求当
\(x=\omega_n^k\) 的时候 \(f(x)\) 的值。\(FFT\) 的分治思想体现在将多项式分为奇次项和偶次项处理。
对于一个 \(n-1\) 次的多项式：
\[f(x) = a_0 + a_1\cdot x + a_2\cdot x^2+...+a_{n-2}\cdot x^{n-2}+a_{n-1}\cdot x^{n-1} \]按照次数的奇偶来分成两组，然后右边提出来一个 \(x\)：
\[\begin{aligned} f(x) &= (a_0+a_2\cdot x^2+a_4\cdot x^4+...+a_{n-2}\cdot x^{n-2}) + (a_1\cdot x+a_3\cdot x^3+a_5\cdot x^5+...+a_{n-1}\cdot x^{n-1})\\ &= (a_0+a_2\cdot x^2+a_4\cdot x^4+...+a_{n-2}\cdot x^{n-2}) + x\cdot (a_1+a_3\cdot x^2+a_5\cdot x^4+a_{n-1}x^{n-2}) \end{aligned} \]分别用奇偶次次项数建立新的函数：
\[\begin{aligned} G(x) &= a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1}\\ H(x) &= a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1} \end{aligned} \]那么原来的 \(f(x)\) 用新函数表示为：
\[f(x)=G\left(x^2\right) + x \times H\left(x^2\right) \]利用偶数次单位根的性质：
\(\omega^i_n = -\omega^{i + \frac{n}{2}}_n\)
\(G\left(x^2\right)\) 和 \(H\left(x^2\right)\) 是偶函数
我们知道在复平面上 \(\omega^i_n\) 和 \(\omega^{i+\frac{n}{2}}_n\) 的 \(G(x^2)\) 的 \(H(x^2)\) 对应的值相同。
得到：
和：
\[\begin{aligned} f(\omega_n^{k+\frac{n}{2}}) &= G(\omega_n^{2k+n}) + \omega_n^{k+\frac{n}{2}} \times H(\omega_n^{2k+n}) \\ &= G(\omega_n^{2k}) - \omega_n^k \times H(\omega_n^{2k}) \\ &= G(\omega_{\frac{n}{2}}^k) - \omega_n^k \times H(\omega_{\frac{n}{2}}^k) \end{aligned} \]因此我们求出了 \(G(\omega_{\frac{n}{2}}^k)\) 和 \(H(\omega_{\frac{n}{2}}^k)\) 后，就可以同时求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\frac{n}{2}})\)。于是对 \(G\) 和 \(H\) 分别递归 \(DFT\) 即可。
考虑到分治 \(DFT\) 能处理的多项式长度只能是 \(2^m(m \in \mathbf{N}^ \ast )\)，否则在分治的时候左右不一样长，右边就取不到系数了。所以要在第一次 \(DFT\) 之前就把序列向上补成长度为 \(2^m(m \in \mathbf{N}^\ast )\)（高次系数补 \(0\)）、最高项次数为 \(2^m-1\) 的多项式。
位逆序置换
原序列的每个数用二进制表示，然后把二进制翻转对称一下，就是最终位置的下标。
可以在 \(O(n\log n)\) 或 \(O(n)\) 的时间内求出每个数变换后的结果
蝶形运算优化
已知 \(G(\omega_{\frac{n}{2}}^k)\) 和 \(H(\omega_{\frac{n}{2}}^k)\) 后，需要使用下面两个式子求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\frac{n}{2}})\)：
\[\begin{aligned} f(\omega_n^k) & = G(\omega_{\frac{n}{2}}^k) + \omega_n^k \times H(\omega_{\frac{n}{2}}^k) \\ f(\omega_n^{k+\frac{n}{2}}) & = G(\omega_{\frac{n}{2}}^k) - \omega_n^k \times H(\omega_{\frac{n}{2}}^k) \end{aligned} \]使用位逆序置换后，对于给定的 \(n\), \(k\)：
\(G(\omega_{\frac{n}{2}}^k)\) 的值存储在数组下标为 \(k\) 的位置，\(H(\omega_{\frac{n}{2}}^k)\) 的值存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。
\(f(\omega_n^k)\) 的值将存储在数组下标为 \(k\) 的位置，\(f(\omega_n^{k+\frac{n}{2}})\) 的值将存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。
因此可以直接在数组下标为 \(k\) 和
\(k + \frac{n}{2}\) 的位置进行覆写，而不用开额外的数组保存值。
再详细说明一下如何借助蝶形运算完成所有段长度为 \(\frac{n}{2}\) 的合并操作：
令段长度为 \(s = \frac{n}{2}\)；
同时枚举序列 \(\{G(\omega_{\frac{n}{2}}^k)\}\) 的左端点 \(l_g = 0, 2s, 4s, \cdots, N-2s\) 和序列 \(\{H(\omega_{\frac{n}{2}}^k)\}\) 的左端点 \(l_h = s, 3s, 5s, \cdots, N-s\)；
合并两个段时，枚举 \(k = 0, 1, 2, \cdots, s-1\)，此时\(G(\omega_\frac{n}{2}^k)\) 存储在数组下标为 \(l_g + k\) 的位置，\(H(\omega_\frac{n}{2}^k)\) 存储在数组下标为 \(l_h + k\) 的位置；
使用蝶形运算求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\frac{n}{2}})\)，然后直接在原位置覆写。
快速傅里叶逆变换
设 \(b_i=\omega_n^{-i}\)，则多项式 \(A\) 在 \(x=b_0,b_1,\cdots,b_{n-1}\) 处的点值表示法为 \(\left\{ A(b_0),A(b_1),\cdots,A(b_{n-1}) \right\}\)。
对 \(A(x)\) 的定义式做一下变换，可以将 \(A(b_k)\) 表示为
\[\begin{aligned} A(b_k)&=\sum_{i=0}^{n-1}f(\omega_n^i)\omega_n^{-ik}=\sum_{i=0}^{n-1}\omega_n^{-ik}\sum_{j=0}^{n-1}a_j(\omega_n^i)^{j}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j-k)}=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i\\ \end{aligned} \]记
\[S\left(\omega_n^a\right)=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i \]当 \(a=0 \pmod{n}\) 时，\(S\left(\omega_n^a\right)=n\)。
当 \(a\neq 0 \pmod{n}\) 时，我们错位相减
也就是说
\[S(\omega_n^a)= \begin{cases} n,a=0\\ 0,a\neq 0 \end{cases} \]那么代回原式
\[A(b_k)=\sum_{j=0}^{n-1}a_jS\left(\omega_n^{j-k}\right)=a_k\cdot n \]也就是说给定点 \(b_i=\omega_n^{-i}\)，则 \(A\) 的点值表示法为
\[\begin{aligned} &\left\{ (b_0,A(b_0)),(b_1,A(b_1)),\cdots,(b_{n-1},A(b_{n-1})) \right\}\\ =&\left\{ (b_0,a_0\cdot n),(b_1,a_1\cdot n),\cdots,(b_{n-1},a_{n-1}\cdot n) \right\} \end{aligned} \]综上所述，我们取单位根为其倒数，对 \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\) 跑一遍 \(FFT\)，然后除以 \(n\) 即可得到 \(f(x)\) 的系数表示。
递归写法
namespace FastFastTle{
il void FFT(qwq *q,int lim,int k){
if(lim==1) return;
for(ri int i=0;i<lim;++i) o[i]=q[i];
ri int m=lim>>1;
for(ri int i=0;i<lim;i++){
if(i&1) q[m+(i>>1)]=o[i];
else q[i>>1]=o[i];
}
FFT(q,m,k),FFT(q+m,m,k);
wn={cos(Pi*2/lim),k*sin(Pi*2/lim)},w={1,0};
for(ri int i=0;i<m;++i,w=w*wn){
x=q[i],y=w*q[i+m];
q[i]=x+y,q[i+m]=x-y;
}
return;
}
}using namespace FastFastTle;
AC·code
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;
namespace Q{
il int rd(){
ri int x=0;ri bool f=0;ri char c=getchar();
while(!isdigit(c)) f|=(c==45),c=getchar();
while(isdigit(c)) x=x*10+(c^48),c=getchar();
return f?-x:x;
}
il void wt(int x){
if(x<0) x=-x,putchar(45);
if(x>=10) wt(x/10);
return putchar(x%10|48),void();
}
} using namespace Q;
cs double Pi=3.141592653589793;
cs int N=2097152;
namespace FastFastTle{
struct qwq{
double r,i;
qwq operator +(cs qwq o)cs{
return (qwq){r+o.r,i+o.i};
}
qwq operator -(cs qwq o)cs{
return (qwq){r-o.r,i-o.i};
}
qwq operator *(cs qwq o)cs{
return (qwq){r*o.r-i*o.i,o.r*i+r*o.i};
}
}wn,w,f[N],g[N],o[N],x,y;
il void FFT(qwq *q,int lim,int k){
if(lim==1) return;
for(ri int i=0;i<lim;++i) o[i]=q[i];
ri int m=lim>>1;
for(ri int i=0;i<lim;i++){
if(i&1) q[m+(i>>1)]=o[i];
else q[i>>1]=o[i];
}
FFT(q,m,k),FFT(q+m,m,k);
wn={cos(Pi*2/lim),k*sin(Pi*2/lim)},w={1,0};
for(ri int i=0;i<m;++i,w=w*wn){
x=q[i],y=w*q[i+m];
q[i]=x+y,q[i+m]=x-y;
}
return;
}
} using namespace FastFastTle;
signed main(){
int n=rd(),m=rd(),lim=1;
while(lim<=n+m) lim<<=1;
for(ri int i=0;i<=n;++i) f[i].r=rd();
for(ri int i=0;i<=m;++i) g[i].r=rd();
FFT(f,lim,1),FFT(g,lim,1);
for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i];
FFT(f,lim,-1);
for(ri int i=0,s;i<=n+m;++i){
wt(f[i].r/lim+0.5),putchar(32);
}
return 0;
}
迭代实现
namespace FastFastTle{
il void FFT(qwq q[],int lim,int k){
for(ri int i=0;i<lim;++i){
if(i<r[i]) swap(q[i],q[r[i]]);
}
for(ri int mid=1;mid<lim;mid<<=1){
wn={cos(Pi/mid),k*sin(Pi/mid)};
for(ri int R=mid<<1,j=0;j<lim;j+=R){
w={1,0};
for(ri int t=0;t<mid;++t,w=w*wn){
x=q[j+t],y=w*q[j+mid+t];
q[j+t]=x+y,q[j+mid+t]=x-y;
}
}
}
return;
}
} using namespace FastFastTle;
AC·code
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;
namespace Q{
il int rd(){
ri int x=0;ri bool f=0;ri char c=getchar();
while(!isdigit(c)) f|=(c==45),c=getchar();
while(isdigit(c)) x=x*10+(c^48),c=getchar();
return f?-x:x;
}
il void wt(int x){
if(x<0) x=-x,putchar(45);
if(x>=10) wt(x/10);
return putchar(x%10|48),void();
}
} using namespace Q;
cs double Pi=3.141592653589793;
cs int N=2097152;
namespace FastFastTle{
struct qwq{
double r,i;
qwq operator +(cs qwq o)cs{
return (qwq){r+o.r,i+o.i};
}
qwq operator -(cs qwq o)cs{
return (qwq){r-o.r,i-o.i};
}
qwq operator *(cs qwq o)cs{
return (qwq){r*o.r-i*o.i,o.r*i+r*o.i};
}
}wn,w,f[N],g[N],x,y;
int l,r[N];
il void FFT(qwq q[],int lim,int k){
for(ri int i=0;i<lim;++i){
if(i<r[i]) swap(q[i],q[r[i]]);
}
for(ri int mid=1;mid<lim;mid<<=1){
wn={cos(Pi/mid),k*sin(Pi/mid)};
for(ri int R=mid<<1,j=0;j<lim;j+=R){
w={1,0};
for(ri int t=0;t<mid;++t,w=w*wn){
x=q[j+t],y=w*q[j+mid+t];
q[j+t]=x+y,q[j+mid+t]=x-y;
}
}
}
return;
}
} using namespace FastFastTle;
signed main(){
int n=rd(),m=rd(),lim=1;
while(lim<=n+m) lim<<=1,l++;
for(ri int i=0;i<=n;++i) f[i].r=rd();
for(ri int i=0;i<=m;++i) g[i].r=rd();
for(ri int i=0;i<lim;++i){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
FFT(f,lim,1),FFT(g,lim,1);
for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i];
FFT(f,lim,-1);
for(ri int i=0,s;i<=n+m;++i){
wt(f[i].r/lim+0.5),putchar(32);
}
return 0;
}